Автор: В. А. Костин
2021 год
$\newcommand{\ctg}{\mathop{\mathrm{ctg}}\nolimits}$
$\newcommand{\tg}{\mathop{\mathrm{tg}}\nolimits}$
$\newcommand{\arctg}{\mathop{\mathrm{arctg}}\nolimits}$
$\newcommand{\degree}{^{\circ}}$
$\renewcommand{\Re}{\mathop{\mathrm{Re}}\nolimits}$
print и вывод ячеек¶При настройках блокнота по умолчанию текстовое представление выходного результата последней команды (кроме инструкций присвоения) в ячейке печатается в вывод ячейки (cell output). Настройки может быть изменены установкой значения InteractiveShell.ast_node_interactivity, InteractiveShell определён в модуле IPython.core.interactiveshell. Возможные значения InteractiveShell.ast_node_interactivity: 'all', 'last', 'last_expr' (по умолчанию), 'none', 'last_expr_or_assign'.
y = 5 + 3
2 + 8
5 + 2
7
Точка с запятой подавляет печать выходного результата последней инструкции.
6;
В Python 3 функция print в Python 3 по умолчанию выводит в стандартный вывод объекты, являющиеся её начальными неименованными аргументами. Стандартный вывод интерпретатора для инструкций в некоторой ячейке с кодом (code cell) в блокноте Jupyter печатается непосредственно после этой ячейки.
print(1, 'text', 3 + 4)
1 text 7
Форматный вывод может осуществляться с помощью метода format переменных и литералов строкового типа.
print('Одна треть равна {0:.3g}, а одна седьмая — {1:.6e}'.format(1 / 3,
1 / 7))
Одна треть равна 0.333, а одна седьмая — 1.428571e-01
Базовые арифметические операции, в целом обозначаются так же, как и в большинстве языков. Круглые скобки могут быть использованы для группировки действий. В частности, сложение выполняется с помощью оператора +.
5 + 2
7
Вычитание выполняется с помощью оператора -.
5 - 2
3
Умножение выполняется с помощью оператора *.
5 * 2
10
Вещественное деление выполняется с помощью оператора /. Результат вещественного деления — всегда число с плавающей точкой.
5 / 2
2.5
Целочисленное деление выполняется с помощью оператора //. Результат — целое число, не превышающее вещественное частное.
5 // 2
2
Остаток от деления может быть получен с помощью оператора %.
5 % 2
1
Оператор ** позволяет выполнять возведение в степень.
5 ** 2
25
Для всех операторов, перечисленных выше, имеются варианты с присвоением: +=, -=, *=. /=, //=. %=, **=.
y = 5
y += 2
y
7
Имеется встроенная поддержка комплексной арифметики.
(2 + 5j) * (2 - 5j)
(29+0j)
Python — язык с динамической типизацией и развитой системой преобразований типов по умолчанию. Тип переменной может быть определён с помощью встроенной функции type.
x = 5 / 2
type(x)
float
x = 5 // 2
type(x)
int
x = 'text'
type(x)
str
Python — строго типизованный язык.
x = 'text'
y = 'something'
try:
y = x + 5
except TypeError:
print('Нельзя прибавить число к строке!')
y
Нельзя прибавить число к строке!
'something'
Python предлагает традиционный набор операторов для контроля потока выполнения. Среди них условный оператор if. Альтернативный блок указывается после оператора else. В случае нескольких условий и альтернатив может быть использован оператор elif. Имеется традиционный набор операторов сравнения ==, <, >, <=, >=, !=.
x = 2
if x == 3:
print('x = 3')
elif x == 2:
print('x = 2')
elif x == 1:
print('x = 1')
else:
print('x is neither 1, 2, or 3')
x = 2
Циклы могут быть организованы с помощью операторов while и for. Внутри цикла могут быть использованы операторы break (для прерывания цикла) и continue (для перехода к следующей итерации).
x = 10
y = x
factorial = 1
while x > 0:
factorial *= x
x -= 1
print('Факториал {0} равен {1}.'.format(y, factorial))
Факториал 10 равен 3628800.
Оператор цикла for используется в паре с оператором in, после которого указывается итерируемый объект (iterable). В качестве такого объекта могут выступать списки, кортежи и словари, а также некоторые другие объекты. Для итерации по арифметической прогрессии может быть использован итератор, конструируемый функцией range. Если range используется с одним целочисленным аргументом, например, range(stop), то итерация проводится от нуля последовательно (через единицу) до stop - 1. Если аргументов 2, range(start, stop), то итерация проводится от start до stop - 1. Третий аргумент может быть использован, для того чтобы указать шаг арифметической прогрессии.
a = 5
b = 15
sum = 0
for i in range(a, b + 1):
sum += i
print('Сумма {0} + {1} + ... + {2} равна {3}.'.format(a, a + 1, b, sum))
Сумма 5 + 6 + ... + 15 равна 110.
Функции могут быть определены с помощью оператора def. Возвращаемый результат может быть указан с помощью оператора return. Функция может возращать несколько значений (то есть возращать кортеж значений).
def f1(x):
return x + 1, 2 * x
f1(3)
(4, 6)
Короткие функции могут быть определены с помощью оператора lambda.
f2 = lambda x: (x + 1, 2 * x)
f2(3)
(4, 6)
Python обладает несколькими встроенными типа составных объектов. В частности, такими типами являются списки и кортежи (неизменяемые списки). Списки конструируются из элементов с помощью квадратных скобок [], кортежи — с помощью круглых ().
x = [1, 2, 5]
type(x)
list
x = (1, 2, 5)
type(x)
tuple
Удобным методом создания списков и кортежей является списковое включение (list comprehension), осуществляемое с помощью операторов квадратных скобок и операторов for и in.
x = [1, 2, 5]
y = [i + 1 for i in x]
y
[2, 3, 6]
Обращение к элементу осуществляется посредством оператора квадратных скобок [].
x = [1, 2, 5]
y = (1, 2, 5)
x[2], y[2]
(5, 5)
Списки и кортежи в Python поддерживают срезы (slicing), то есть обращение к подсписку или подкортежу, заданному арифметической прогрессией индексов.
x = [i for i in range(10)]
print('Исходный список: {0}'.format(x))
print('Его слайс (срез) [3:9:2]: {0}'.format(x[3:9:2]))
Исходный список: [0, 1, 2, 3, 4, 5, 6, 7, 8, 9] Его слайс (срез) [3:9:2]: [3, 5, 7]
x = [i for i in range(10)]
print('Исходный список: {0}'.format(x))
print('Список в противоположном порядке: {0}'.format(x[::-1]))
Исходный список: [0, 1, 2, 3, 4, 5, 6, 7, 8, 9] Список в противоположном порядке: [9, 8, 7, 6, 5, 4, 3, 2, 1, 0]
Любой объект в Python может быть изменяемым (mutable) или неизменяемым (immutable). Числа и строки — примеры неизменяемых объектов. Неизменяемыми являются объекты числовых типов, булевы значения, строки, неизменяемые последовательности байтов (тип bytes), кортежи (тип tuple), неизменяемые (замороженные) множества (тип frozenset). Остальные объекты, например, списки, являются изменяемыми объектами. Присваивание переменной изменяемого или неизменяемого объекта носит различный характер.
# Числа — неизменяемые (immutable) объекты
x = 4 # Имя «x» связывается с неизменяемым объектом «4»
y = x # Имя «y» связывается с неизменяемым объектом «4»
y += 1 # Имя «y» связывается с неизменяемым объектом «5», значение x не
# меняется
x
4
# Списки — изменяемые (mutable) объекты
# объектами
x = [4] # Выделяется память, в которую записывается список «[4]», имя «x»
# связывается с адресом списка в памяти
y = x # Имя «y» связывается с тем же адресом, копии объекта не создаётся
y[0] += 1 # y[0] и x[0] — синонимы
x
[5]
При работе с Python в неинтерактивном режиме для библиотеки Numpy рекомедуется использовать стандартную инструкцию импорта import numpy или для сокращения import numpy as np и обращаться к фукнциям и классам библиотеки через модуль, например, numpy.arange или np.arange. При интерактивной работе, например, в Jupyter, для подключения объектов Numpy и Matplotlib удобно воспользоваться модулем pylab и импортировать все объекты в глобальное пространство имён. При этом функции и классы оказываются доступны без указания библиотеки, то есть используется arange вместо numpy.arange.
from pylab import *
Библиотека вводит новый тип составного объекта — numpy.ndarray. Этот тип позволяет хранить и проводить операции с многомерными массивами однородных (то есть однотипных) данных. Функция array позволяет создать такой массив из других составных объектов, например, из простых и вложенных списков или кортежей.
x = array([[0, 1, 2], [3, 4, 5]])
x
array([[0, 1, 2],
[3, 4, 5]])
type(x)
numpy.ndarray
Тип элемента массива содержится в поле dtype.
x.dtype
dtype('int32')
Поле ndim содержит число измерений массива.
x.ndim
2
Поле size содержит число элементов массива.
x.size
6
Поле shape содержит кортеж с размерами массива вдоль различных измерений
x.shape
(2, 3)
Доступ к отдельному элементу может быть осуществлён посредством оператора квадратных скобок [].
x[1, 2], x[1][2]
(5, 5)
Массив содержащий арифметическую прогрессию может быть создан с помощью функции arange. Смысл первых трёх аргументов функции аналогичен смыслу аргументов функции range, но в качестве аргументов могут выступать не только целые числа, но и числа с плавающей точкой.
arange(1., 4., 0.3)
array([1. , 1.3, 1.6, 1.9, 2.2, 2.5, 2.8, 3.1, 3.4, 3.7])
Функция linspace также как функция arange выдаёт массив, содержащий арифметическую прогрессию, но принимает в качестве аргументов начальное значение, конечное значение и размер выходного массива.
linspace(2, 4, 5)
array([2. , 2.5, 3. , 3.5, 4. ])
Функция meshgrid конструирует координатные матрицы из координатных векторов.
x = arange(1, 4)
y = arange(1, 4)
X, Y = meshgrid(x, y)
print('Матрица X:\n{0}'.format(X))
print('Матрица Y:\n{0}'.format(Y))
Матрица X: [[1 2 3] [1 2 3] [1 2 3]] Матрица Y: [[1 1 1] [2 2 2] [3 3 3]]
Функция zeros создаёт массив заданных формы и типа, заполненный нулями.
zeros((2, 3), dtype=double)
array([[0., 0., 0.],
[0., 0., 0.]])
Функция full создаёт массив с одинаковыми значениями всех элементов.
full((2, 3), 5.)
array([[5., 5., 5.],
[5., 5., 5.]])
Функция empty создаёт пустой массив (то есть при его создания память не перезаписывается какими-либо заданными значениями).
empty((2, 3), dtype=double)
array([[5., 5., 5.],
[5., 5., 5.]])
Функции zeros_like, full_like и empty_like действуют как функции zeros, full и empty, но принимают в качестве аргумента вместо формы другой массив у которого копирует форму и тип (если тот не указывается явно среди аргументов).
x = zeros((2, 3), dtype=int)
empty_like(x)
array([[475699600, 477, 0],
[ 0, 1, 32767]])
Так как массивы Numpy являются изменяемыми объектами, присваивание не создаёт копии массива. Для создания копии массива используется функция copy.
x = arange(10)
y = x
y[5] = -100
x
array([ 0, 1, 2, 3, 4, -100, 6, 7, 8, 9])
x = arange(10)
y = copy(x)
y[5] = -100
x
array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
Для массивов Numpy переопределены арифметические операции и элементарные математические функции. В большистве случаев операции и функции применяются поэлементно. При этом для бинарных операций оба операнда должны иметь одинаковую форму либо один из операндов должен быть скаляром. При выполнении операций тип элемента массива может меняться (например, при делении целочисленных массивов чисел или вычислении трансцендентых функций от таких массивов).
x = array([[0, 1, 2], [3, 4, 5]])
x
array([[0, 1, 2],
[3, 4, 5]])
y = x + 1
y
array([[1, 2, 3],
[4, 5, 6]])
x * y
array([[ 0, 2, 6],
[12, 20, 30]])
z = sin(x)
z
array([[ 0. , 0.84147098, 0.90929743],
[ 0.14112001, -0.7568025 , -0.95892427]])
x.dtype, z.dtype
(dtype('int32'), dtype('float64'))
Массивы Numpy поддерживают срезы (слайсинг) для операций с подмассивами.
x = arange(10)
x
array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
Для положительных c срез a:b:c выбирает каждый c-й элемент с номерами ниже b, начиная с элемента номером a, при этом значения a и b берутся по модулю размера массива в соответствующем измерении.
x[3:9:2]
array([3, 5, 7])
Для отрицательных c срез a:b:c выбирает каждый c-й элемент в обратном порядке с номерами выше b, начиная с элемента номером a.
x[9:3:-2]
array([9, 7, 5])
Каждое из числе a, b, c в срезе a:b:c может быть опущено. Значение a и c по умолчанию равны 0 и 1 соответственно, значение b по умолчанию равно размеру массива по соответствующему измерению. Двоеточие в конце среза также может опускаться.
x[::-1]
array([9, 8, 7, 6, 5, 4, 3, 2, 1, 0])
Одиночное двоеточие обозначает все элементы в каком-то измерении. Его удобно использовать, чтобы выбирать проекции многомерного массива, например, выбирать некоторые столбцы в двумерном массиве (матрице).
y = zeros((5, 7))
y[::2, ::2] = 1
y
array([[1., 0., 1., 0., 1., 0., 1.],
[0., 0., 0., 0., 0., 0., 0.],
[1., 0., 1., 0., 1., 0., 1.],
[0., 0., 0., 0., 0., 0., 0.],
[1., 0., 1., 0., 1., 0., 1.]])
y[:, 2]
array([1., 0., 1., 0., 1.])
При аналогичном выборе строк в матрице запятую с двоеточием в конце можно опустить.
y[2]
array([1., 0., 1., 0., 1., 0., 1.])
Библиотека Matplotlib предлагает два подхода к построению рисунков императивный (схожий с Matlab) и объектно-ориентированный. Для простой графики, по всей видимости, более удобен императивный подход, реализуемый с помощью функций модуля matplotlib.pyplot. При работе с Python в неинтерактивном режиме для его подключения рекомедуется использовать стандартную инструкцию импорта import matplotlib.pyplot as plt и обращаться к фукнциям и классам библиотеки через модуль, например, plt.plot. При интерактивной работе, например, в Jupyter, для подключения объектов Numpy и Matplotlib удобно воспользоваться модулем pylab и импортировать все объекты в глобальное пространство имён. При этом функции и классы оказываются доступны без указания библиотеки, то есть используется plot вместо plt.plot.
from pylab import *
Так называемая «магическая» команда %matplotlib позволяет выбрать тип вывода (бэкенд) рисунков. Параметр inline указывает на то, что рисунки Matplotlib должны быть встроены как растровые изображения, для интерактивных рисунков на JavaScript можно использовать параметр notebook вместо inline. Команда %matplotlib --list выводит список возможных типов вывода (бэкендов).
%matplotlib inline
Переменная rcParams содержит различные настройки рисунков. Доступ к настройкам осуществляется по ключам. Например, ключи figure.figsize и figure.dpi соответствуют размеру рисунков по умолчанию в дюймах и разрешению рисунков в точках на дюйм.
rcParams['figure.figsize'] = (9.6, 7.2)
rcParams['figure.dpi'] = 100
Один из наиболее простых способов построить график функции или параметрически заданную кривую в Matplotlib — использовать функцию plot. Функция принимает в качестве аргумента массивы абсцисс и ординат точек кривой. Рисунок можно снабдить заголовком с помощью функции title, а подписать оси можно с помощью функций xlabel и ylabel. Добавить подпись в произвольном месте можно с помощью команды text. Подписи могут содержать формулы в формате, схожем с LaTeX.
x = arange(0., 10., 0.1)
y = sin(x)
plot(x, y)
xlabel(r'$x$') # Подписи по осям могут содержать формулы
ylabel(r'$\sin x$') # r перед строкой позволяет не экранировать специальные
# символы, в частности, обратную косую черту \
text(3., 0.5, r'Кривая 1')
title('График синуса'); # Точка с запятой подавляет вывод последней команды
Функция plot также может принимать аргументы, контролирующие параметры отображения кривой.
plot(x, y, linewidth=4.0)
xlabel(r'$x$')
ylabel(r'$\sin x$')
title('График синуса');
Третьим аргументов функция plot может принимать стилевую строку, указывающую цвет кривой, её тип и вид маркеров.
plot(x, y, 'ro')
xlabel(r'$x$')
ylabel(r'$\sin x$')
title('График синуса');
Функция plot может несколько наборов абсцисс, ординат и стилевых строк для построения нескольких графиков в одних осях.
plot(x, y, x[::9], y[::9], 'ro', x, cos(x))
xlabel(r'$x$')
ylabel(r'$\sin x$ и $\cos x$')
title('Графики синуса и косинуса');
Несколько графиков можно построить, вызвав функцию plot несколько раз.
plot(x, y)
plot(x[::9], y[::9], 'ro')
plot(x, cos(x))
xlabel(r'$x$')
ylabel(r'$\sin x$ и $\cos x$')
title('Графики синуса и косинуса');
Функция fill аналогична функции plot, но вместо графика или точек рисует внутреннюю область графика, заполненную цветом. В примере ниже инструкция axis('equal') устанавливает одинаковые масштабы по горизонтальной и вертикальной осям.
t = linspace(0., 2 * pi, 201)
fill(cos(t), sin(t), color='pink')
fill([2, 2, 4, 4], [-1, 1, 1, -1], color='purple')
axis('equal')
xlabel(r'$x$')
xlabel(r'$y$')
title('Круг и квадрат');
Несколько координатных плоскостей (осей, осевых коробок, axes) на одном рисунке можно получить с помощью функции subplot. Команда subplot(nrows, ncols, index) создаёт виртуальную сетку из координатных плоскостей с nrows строками и ncols столбцами и возвращает координатную плоскость с индексом index из этой сетки. Общий заголовок для всех координатных плоскостей можно получить с помощью функции suptitle. Функция tight_layout может использовать для корректировки размеров осей тем, чтобы подписи по осям не перекрывали других подрисунков и не выходили за пределы всего рисунка. Аргумент rect указывает на положение граничной рамки, содержащей все подрисунки, но не главный заголовок.
x = arange(0., 10., 0.1)
subplot(2, 1, 1) # Первая панель (первый «подрисунок»)
plot(x, sin(x))
xlabel(r'$x$')
ylabel(r'$\sin x$')
subplot(2, 1, 2) # Вторая панель (второй «подрисунок»)
plot(x, cos(x))
xlabel(r'$x$')
ylabel(r'$\cos x$')
suptitle(r'Графики синуса и косинуса')
tight_layout(rect=[0, 0, 1, 0.95]);
Функцию streamplot удобно использовать для построения фазовых портретов, силовых линий векторных полей или линий тока на плоскости. В простейшем варианте функция принимает четыре аргумента: первые два задают координатные матрицы, два следующих задают две компоненты векторного поля. Одним из возможных дополнительных параметров является density, определяющий характерную плотность отображаемых линий (то есть обратное расстояние между линиями). Далее приводится пример построения фазовой портрета с использованием функции streamplot для системы $\ddot\varphi + 0.3\dot\varphi + \sin\varphi = 0$.
# Создание сеток переменных 201 на 201 элемента
phi, dphi = meshgrid(linspace(-pi, pi, 201), linspace(-3., 3., 201))
# Вычисление вертикальной компоненты векторного поля (горизонтальная
# компонента равна dphi и её не нужно специально вычислять)
ddphi = -sin(phi) - 0.3 * dphi
# Построение фазовой плоскости: первые два аргумента задают сетку значений
# абсцисс и ординат, следующие два аргумента задают горизонтальную и
# вертикальную компоненты векторного поля
streamplot(phi, dphi, dphi, ddphi, density=1.5, color='k')
xlabel(r'$\varphi$')
ylabel(r'$\dot\varphi$')
title(r'Фазовый портрет системы $\ddot\varphi + 0.3\dot\varphi + ' +
r'\sin\,\varphi = 0$');
Линии уровня некоторой функции двух переменных можно построить с помощью функции contour. В одном из вариантов функция принимает первыми аргументами две координатные матрицы и матрицу значений функции. Четвёртый аргумент может задать количество линий уровня или упорядоченный по возрастанию набор значений фукнции для которых строятся линии уровня. Функция clabel может быть использована для расстановки меток значений уровня на линиях. В этой функции параметр inline контролирует, стираются ли линии уровня под метками, параметр fontsize управляет размером шрифта меток. Рисунок также можно снабдить шкалой с помощью функции colorbar.
x = linspace(-2., 2., 200)
y = linspace(-2., 2., 200)
X, Y = meshgrid(x, y)
Z = X * exp(-X ** 2 - 1.5 * Y ** 2) + 0.2 * Y
v = linspace(-0.5, 0.5, 21)
contours = contour(X, Y, Z, v)
gca().set_aspect('equal') # Установка одинаковых масштабов по вертикальной
# и горизонтальной осям (здесь необязательно)
clabel(contours, inline=True, fontsize=8)
colorbar(contours, orientation='vertical', label=r'$z$')
xlabel(r'$x$')
ylabel(r'$y$')
title(r'Линии уровня $z = e^{-x^2 - 1.5y^2} + 0.2y$');
Функция contourf работает аналогично contour но строит не сами линии уровня, а закрашивает разными цветами области между линиями уровня. Такой рисунок также можно снабдить шкалой с помощью функции colorbar.
x = linspace(-2., 2., 200)
y = linspace(-2., 2., 200)
X, Y = meshgrid(x, y)
Z = X * exp(-X ** 2 - 1.5 * Y ** 2) + 0.2 * Y
v = linspace(-0.5, 0.5, 21)
contours = contourf(X, Y, Z, v)
gca().set_aspect('equal') # Установка одинаковых масштабов по вертикальной
# и горизонтальной осям (здесь необязательно)
colorbar(contours, orientation='vertical', label=r'$z$')
xlabel(r'$x$')
ylabel(r'$y$')
title(r'Тепловая карта $z = e^{-x^2 - 1.5y^2} + 0.2y$');
Упомянутой выше в разделе 5.3.5 функции contourf может использоваться для построения векторной тепловой карты, представляющей набор областей равных цветов. Для построения растровой тепловой карты можно применять функцию imshow. В этой функции параметр origin указывает, верхнему или нижнему левому углу панели соответствует первый (с индексами 0, 0) элемент изображаемого массива. Обычно этот параметр должен быть установлен в 'lower' при построении тепловых карт функций, заданных матрицами в традиционных декартовых координатах. Параметр extent указывает координаты левого нижнего и правого верхнего углов тепловой карты.
xlims = [-2., 2.]
ylims = [-2., 2.]
x = linspace(xlims[0], xlims[1], 200)
y = linspace(ylims[0], ylims[1], 200)
# Переход к средним точкам, чтобы значение функции соответствовало центру
# пикселя. При высоком разрешении (большом количестве точек) необязательно
x = convolve(x, [0.5, 0.5], 'valid')
y = convolve(y, [0.5, 0.5], 'valid')
X, Y = meshgrid(x, y)
Z = X * exp(-X ** 2 - 1.5 * Y ** 2) + 0.2 * Y
im = imshow(Z, interpolation='bilinear', origin='lower', aspect='auto',
extent=(xlims[0], xlims[1], ylims[0], ylims[1]))
gca().set_aspect('equal') # Установка одинаковых масштабов по вертикальной
# и горизонтальной осям (здесь необязательно)
colorbar(im, orientation='vertical', label=r'$z$')
xlabel(r'$x$')
ylabel(r'$y$')
title(r'Тепловая карта $z = e^{-x^2 - 1.5y^2} + 0.2y$');
control_theory¶Plant¶class Plant(seed_value=0)
Класс Plant — является базовым интерфейсным классом для представления объекта управления, задаваемого с помощью системы обыкновенных дифференциальных уравнений
$$\dot{\mathbf x} = \mathbf f(\mathbf x, \mathbf u, t),$$
где $\mathbf x$ — векторная переменная, определяющеая состояние объекта (элемент пространства состояний), $\mathbf u$ — векторное управление, функция $\mathbf f$ полностью определяет модель объекта управления.
Экземпляр базового класса Plant представляет пустой (тривиальный) объект управления без переменных, параметров и управления. Экземпляры классов-наследников представляют специфичные нетривиальные объекты управления.
Принимаемые аргументы
| Название | Тип | Описание | Значение по умолчанию |
|---|---|---|---|
seed_value |
int |
Зерно генератора случайных чисел для инициализации параметров объекта управления (если такие параметры есть). | 0 |
Plant.PLANT_DIM, Plant.CONTROL_DIM, Plant.OTYPE¶Класс и его наследники содержит следующие поля класса (то есть поля, общие для всех экземпляров класса).
| Название | Тип | Описание |
|---|---|---|
PLANT_DIM |
int |
Размерность состояния объекта (то есть вектора $\mathbf x$). |
CONTROL_DIM |
int |
Размерность управления (то есть вектора $\mathbf u$). |
OTYPE |
str |
Название (строковый идентификатор) типа объекта управления. |
Plant.__call__¶f = Plant.__call__(x, u=zeros(CONTROL_DIM), t=0.)
Экзмепляр класса является функтором и при вызове вычисляет значение временной производной $\dot {\mathbf x}$, иными словами вычисляет функцию $\mathbf f$.
Принимаемые аргументы
| Название | Тип | Описание | Значение по умолчанию |
|---|---|---|---|
x |
array_like, shape=(PLANT_DIM,) |
Вектор состояния объекта $\mathbf x$. | |
u |
array_like, shape=(CONTROL_DIM,) |
Вектор управления $\mathbf u$ (по умолчанию — нулевой). | zeros(CONTROL_DIM) |
t |
number |
Значение времени $t$. | 0. |
Возвращаемые значения
| Название | Тип | Описание |
|---|---|---|
f |
ndarray, shape=(PLANT_DIM,), dtype='double' |
Вектор временной производной $\mathbf f(\mathbf x, \mathbf u, t)$ состояния объекта. |
Plant.__str__¶plant_str = Plant.__str__()
Для экзмепляра класса определена функция преобразования в строку __str__.
Возвращаемые значения
| Название | Тип | Описание |
|---|---|---|
plant_str |
str |
Строковое представление экземпляра класса. |
Pendulum¶Класс Pendulum — потомок класса Plant для представления маятника на каретке, описываемого системой $(8)$.
Pendulum.omega, Pendulum.alpha, Pendulum.lambd, Pendulum.mu, Pendulum.nu¶Экземпляр класса содержит следующие поля, соответствующие параметрам системы $(8)$. При инициализации объекта связываются со случайными действительными числами.
| Название | Тип | Описание |
|---|---|---|
omega |
float |
Параметр $\Omega$, значение при инициализации равномерно распределено в $[0.3, 1.8)$. |
alpha |
float |
Параметр $\alpha$, значение при инициализации равномерно распределено в $[0.3, 0.8)$. |
lambd |
float |
Параметр $\lambda$, значение lambd - nu при инициализации равномерно распределено в $[0.01, 0.11)$. |
mu |
float |
Параметр $\mu$, значение mu - alpha * nu при инициализации равномерно распределено в $[0.01, 0.11)$. |
nu |
float |
Параметр $\nu$, значение при инициализации равномерно распределено в $[0.01, 0.11)$. |
Control¶class Control()
Класс Control — является базовым интерфейсным классом для представления стратегии управления, задаваемой с помощью системы уравнений
$$\mathbf u = \mathbf g(\mathbf x, \mathbf v, t),$$
$$\dot{\mathbf v} = \mathbf h(\mathbf x, \mathbf v, t),$$
где $\mathbf x$ — векторная переменная, определяющеая состояние объекта (элемент пространства состояний), $\mathbf u$ — векторное управление, $\mathbf v$ — внутренние («скрытые») переменные, характеризующие состояние регулятора, функция $\mathbf h$ определяет внутреннюю динамику скрытых переменных регулятора, функция $\mathbf g$ связывает управление со состояниями объекта управления и регулятора.
Экземпляр базового класса Control представляет пустой регулятор без переменных. Экземпляры классов-наследников представляют специфичные нетривиальные регуляторы.
Control.PLANT_DIM, Control.CONTROL_DIM, Control.HIDDEN_DIM, Control.CTYPE¶Класс и его наследники содержит следующие поля класса (то есть поля, общие для всех экзмепляров класса).
| Название | Тип | Описание |
|---|---|---|
PLANT_DIM |
int |
Размерность состояния объекта (то есть вектора $\mathbf x$). |
CONTROL_DIM |
int |
Размерность управления (то есть вектора $\mathbf u$). |
HIDDEN_DIM |
int |
Размерность состояния регулятора (то есть вектора $\mathbf v$). |
CTYPE |
str |
Название (строковый идентификатор) типа регулятора. |
Control.__call__¶g = Control.__call__(x, v, t=0.)
Экзмепляр класса является функтором и при вызове вычисляет значение временной производной $\dot {\mathbf v}$, иными словами вычисляет функцию $\mathbf h$.
Принимаемые аргументы
| Название | Тип | Описание | Значение по умолчанию |
|---|---|---|---|
x |
array_like, shape=(PLANT_DIM,) |
Вектор состояния объекта $\mathbf x$. | |
v |
array_like, shape=(HIDDEN_DIM,) |
Вектор состояния регулятора $\mathbf v$. | |
t |
number |
Значение времени $t$. | 0. |
Возвращаемые значения
| Название | Тип | Описание |
|---|---|---|
g |
ndarray, shape=(HIDDEN_DIM,), dtype='double' |
Вектор временной производной $\mathbf h(\mathbf x, \mathbf v, t)$ состояния регулятора. |
Control.control¶u = Control.control(x, v, t=0.)
Метод control вычисляет управление $\mathbf u$ по состояниям объекта и регулятора, иными словами вычисляет функцию $\mathbf g$.
Принимаемые аргументы
| Название | Тип | Описание | Значение по умолчанию |
|---|---|---|---|
x |
array_like, shape=(PLANT_DIM,) |
Вектор состояния объекта $\mathbf x$. | |
v |
array_like, shape=(HIDDEN_DIM,) |
Вектор состояния регулятора $\mathbf v$. | |
t |
number |
Значение времени $t$. | 0. |
Возвращаемые значения
| Название | Тип | Описание |
|---|---|---|
u |
ndarray, shape=(CONTROL_DIM,), dtype='double' |
Вектор управления $\mathbf u$. |
Control.__str__¶control_str = Control.__str__()
Для экзмепляра класса определена функция преобразования в строку __str__.
Возвращаемые значения
| Название | Тип | Описание |
|---|---|---|
control_str |
str |
Строковое представление экземпляра класса. |
ZeroControl¶class ZeroPendulumControl(plant_dim, control_dim)
Класс ZeroPendulumControl — потомок класса Control для представления нулевого управления для маятником на тележке Pendulum.
Принимаемые аргументы
| Название | Тип | Описание | Значение по умолчанию |
|---|---|---|---|
plant_dim |
int |
Размерность состояния объекта управления. | 1 |
control_dim |
int |
Размерность управления. | 1 |
StateFeedbackControl¶class StateFeedbackControl(plant_dim, control_dim, fun)
Класс StateFeedbackControl — потомок класса ZeroControl для представления управления по состоянию объекта.
Принимаемые аргументы
| Название | Тип | Описание | Значение по умолчанию |
|---|---|---|---|
plant_dim |
int |
Размерность состояния объекта управления. | 1 |
control_dim |
int |
Размерность управления. | 1 |
fun |
u = callable(x); x: arraylike, shape=(plant_dim,); u: arraylike, (control_dim,) или None |
Функция, задающая стратегию управления: принимает текущее состояние объекта x и возращает управление u. Если передаётся None, то функция считается возвращающей нулевой вектор (поведение класса не отличается от ZeroControl). |
None |
StateFeedbackControl.fun¶Экземпляр класса содержит поле с функцией, задающей стратегию управления.
| Название | Тип | Описание |
|---|---|---|
fun |
u = callable(x); x: arraylike, shape=(PLANT_DIM,); u: arraylike, (CONTROL_DIM,) |
Функция, задающая стратегию управления: принимает текущее состояние объекта x и возращает управление u. |
StateTimeFeedbackControl¶class StateTimeFeedbackControl(plant_dim, control_dim, fun)
Класс StateFeedbackControl — потомок класса ZeroControl для представления комбинированного управления по состоянию объекта и моменту времени (программного управления).
Принимаемые аргументы
| Название | Тип | Описание | Значение по умолчанию |
|---|---|---|---|
plant_dim |
int |
Размерность состояния объекта управления. | 1 |
control_dim |
int |
Размерность управления. | 1 |
fun |
u = callable(x, t); x: arraylike, shape=(plant_dim,); t: number; u: arraylike, (control_dim,) или None |
Функция, задающая стратегию управления: принимает текущее состояние объекта x и время t и возращает управление u. Если передаётся None, то функция считается возвращающей нулевой вектор (поведение класса не отличается от ZeroControl). |
None |
StateTimeFeedbackControl.fun¶Экземпляр класса содержит поле с функцией, задающей стратегию управления.
| Название | Тип | Описание |
|---|---|---|
fun |
u = callable(x, t); x: arraylike, shape=(PLANT_DIM,); t: number; u: arraylike, (CONTROL_DIM,) |
Функция, задающая стратегию управления: принимает текущее состояние объекта x и время t и возращает управление u. |
LinearStateControl¶class LinearStateControl(plant_dim, control_dim, k=None)
Класс LinearStateControl — потомок класса ZeroControl для представления линейного управления по состоянию.
Принимаемые аргументы
Принимаемые аргументы
| Название | Тип | Описание | Значение по умолчанию |
|---|---|---|---|
plant_dim |
int |
Размерность состояния объекта управления. | 1 |
control_dim |
int |
Размерность управления. | 1 |
k |
arraylike, shape=(control_dim, plant_dim) или shape=(plant_dim,), если control_dim равно 1, или None |
Матрица или вектор $\mathbf k$ коэффициентов линейного закона управления. Задание None эквивалентно заданию нулевой матрицы с shape=(control_dim, plant_dim). |
None |
LinearStateControl.k¶Экземпляр класса содержит поле с матрицей $\mathbf k$ коэффициентов линейного закона управления.
| Название | Тип | Описание |
|---|---|---|
k |
ndarray, shape=(CONTROL_DIM, PLANT_DIM), dtype='double' |
Матрица коэффциентов $\mathbf k$ линейного закона управления (матрица умножается на вектор слева). |
LinearSinPendulumControl¶class LinearSinPendulumControl(k=zeros(4, dtype='double'))
Класс LinearSinPendulumControl — потомок класса Control для представления управления $u = k_0 \sin(x_0) + k_1 x_1 + k_2 x_2 + k_3 x_3$ для маятника на каретке.
Принимаемые аргументы
| Название | Тип | Описание | Значение по умолчанию |
|---|---|---|---|
k |
array_like, shape=(4,) |
Вектор $\mathbf k$ коэффициентов линейного закона управления. | zeros(4, dtype='double') |
LinearSinPendulumControl.k¶Экземпляр класса содержит поле с вектором $\mathbf k$ коэффициентов линейного закона управления.
| Название | Тип | Описание |
|---|---|---|
k |
ndarray, shape=(4,), dtype='double' |
Вектор коэффциентов $\mathbf k$. |
rk¶tout, y = rk(rhs, y0, t, N=1, stop_condition=None)
Функция rk реализует классический метод Рунге — Кутты четвёртого порядка для решения обыкновенных дифферециальных уравений.
Принимаемые аргументы
| Название | Тип | Описание | Значение по умолчанию |
|---|---|---|---|
rhs |
dydt1 = callable(y1, t1); y1: ndarray, dtype=float or dtype=complex; t1: number; dydt1: ndarray, shape=y1.shape |
Функция, задающая систему обыкновенных уравнений: для массива y1 значений зависимых переменных (координат) и независимой переменной (времени) t1 возвращает массив dydt1 значений производных координат по времени. |
|
y0 |
arraylike, shape=y1.shape |
Массив начальных значений координат. | |
t |
arraylike, ndim=1 |
Сортированный по возрастанию или убыванию массив значений времени, для которых вычисляются значения координат. Первый элемент массива отвечает времени для которого задаётся начальное значение. Если массив сортирован по убыванию, расчёт ведётся в обратном времени (то есть в сторону уменьшения времени). | |
N |
int |
Число шагов метода Рунге — Кутты, приходящееся на один интервал, образуемый соседними элементами из t. |
1 |
stop_condition |
None или stop_flag = callable(y1, t1); y1: ndarray, dtype=float or dtype=complex; t1: number; stop_flag: bool |
Функция, задающая критерий остановки, если для текущих значений y1 и t1 выдаёт True, выполнение останавливается и выдаются предшествующие рассчитанные значения времени и переменного вектора. Если задано None выполнение останавливается, когда пройден весь массив времён t. |
None |
Возращаемые значения
| Название | Тип | Описание |
|---|---|---|
tout |
ndarray, ndim=y1.ndim + 1, ndim=1, shape[0] <= t.shape[0] |
Массив, содержащий значения времён, для которых были рассчитаны значения координат. Если функция stop_condition не задана или не выдавала True в ходе расчёта, выводится копия t. Иначе выводится обрезанный вариант массива t, при этом последний элемент массива отвечает времени, на котором условие остановки было выполнено. |
y |
ndarray, ndim=y1.ndim + 1, shape=(tout.shape[0],) + y1.shape |
Многомерный массив, содержащий рассчитанные значения координат для времён из tout. |
pcrhs¶rhs = pcrhs(p, c)
Функция pcrhs для заданных объекта и регулятора создаёт функцию (функциональный объект), которая рассчитывает правую часть системы обыкновенных дифференциальных уравнений, описывающих совместную эволюцию состояний объекта и регулятора во времени
$$\dot{\mathbf X} = \mathbf F(\mathbf X, t),$$
где
$\mathbf X = \begin{bmatrix} \mathbf x \\ \mathbf v\end{bmatrix}$ — объединённый вектор состояний объекта и регулятора, $\mathbf F(\mathbf X, t) = \begin{bmatrix} \mathbf f[\mathbf x, \mathbf g(\mathbf x, \mathbf v, t), t] \\ \mathbf h(\mathbf x, \mathbf v, t) \end{bmatrix}$. Эта система эквивалентна системе приведённых ранее уравнений
$$\dot{\mathbf x} = \mathbf f(\mathbf x, \mathbf u, t),$$
$$\mathbf u = \mathbf g(\mathbf x, \mathbf v, t),$$
$$\dot{\mathbf v} = \mathbf h(\mathbf x, \mathbf v, t).$$
Принимаемые аргументы
| Название | Тип | Описание | Значение по умолчанию |
|---|---|---|---|
p |
Plant |
Объект управления. | |
c |
Control, PLANT_DIM=p.PLANT_DIM, CONTROL_DIM=p.CONTROL_DIM |
Регулятор. |
Возращаемые значения
| Название | Тип | Описание |
|---|---|---|
rhs |
dydt1 = callable(y1, t1); y1: ndarray, dtype=float or dtype=complex, shape=(p.PLANT_DIM + c.HIDDEN_DIM,); t1: number; dydt1: ndarray, shape=y1.shape |
Функция, вычисляющая зависимость $\mathbf F(\mathbf X, t)$. |
control_output¶ufun = control_output(c)
Функция control для заданного регулятора создаёт функцию (функциональный объект), которая рассчитывает управление по состояниям объекта и управления, $\mathbf u(\mathbf X, t) = \mathbf g(\mathbf x, \mathbf v, t)$.
Принимаемые аргументы
| Название | Тип | Описание | Значение по умолчанию |
|---|---|---|---|
c |
Control, PLANT_DIM=p.PLANT_DIM, CONTROL_DIM=p.CONTROL_DIM |
Регулятор. |
Возращаемые значения
| Название | Тип | Описание |
|---|---|---|
ufun |
u = callable(y1, t1); y1: ndarray, dtype=float or dtype=complex, shape=(c.PLANT_DIM + c.HIDDEN_DIM,); t1: number; u: ndarray, shape=(c.CONTROL_DIM,) |
Функция, вычисляющая зависимость $\mathbf u(\mathbf X, t)$. |
integrate¶t, y = integrate(p, c, x0, v0, dt, T, N=1, method=rk, return_control=False, stop_condition=None)
или
t, y, u = integrate(p, c, x0, v0, dt, T, N=1, method=rk, return_control=False, stop_condition=None)
Функция integrate для заданных объекта и регулятора рассчитывает эволюцию их состояний во времени посредством численного интегрирования соответствующих обыкновенных дифференциальных уравнений.
Принимаемые аргументы
| Название | Тип | Описание | Значение по умолчанию |
|---|---|---|---|
p |
Plant |
Объект управления. | |
c |
Control, PLANT_DIM=p.PLANT_DIM, CONTROL_DIM=p.CONTROL_DIM |
Регулятор. | |
x0 |
array_like, shape=(p.PLANT_DIM,) |
Начальное состояние объекта. | |
v0 |
array_like, shape=(c.HIDDEN_DIM,) |
Начальное состояние регулятора. | |
dt |
number |
Шаг по времени в выходных массивах. Если задаётся отрицательным, то расчёт ведётся в обратном времени. | |
T |
number |
Максимальное (для положительного dt) или минимальное (для отрицательного dt) время, до которого рассчитывается эволюция. |
|
N |
int |
Число шагов метода интегрирования, приходящееся на один шаг по времени в выходных массивах. | 1 |
method |
callable |
Функция, реализующая метод интегрирования с сигнатурой, как у rk. |
rk |
return_control |
Если True, то возвращает массив значений управления помимо массивов времени и состояний. |
False |
|
stop_condition |
None or b1 = callable(y1, t1); y1: ndarray, shape=(p.PLANT_DIM + c.HIDDEN_DIM,); t1: number; b1: bool |
Необязательная функция координат и времени, задающая условие остановки расчёта, если функция возвращает True после шага, расчёт останавливается. |
None |
Возращаемые значения
| Название | Тип | Описание |
|---|---|---|
t |
ndarray, ndim=1 |
Эквидистантный массив значений времени, для которых выводятся состояния объекта и регулятора. Если функция stop_condition не задана или не сработала в ходе расчёта, то форма выходного массива shape=(int(ceil(T/dt)),). Если функция stop_condition сработала, то последний элемент массива отвечает времени, на котором условие остановки было выполнено. |
y |
ndarray, shape=(t.shape[0], p.PLANT_DIM + c.HIDDEN_DIM) |
Двумерный массив, содержащий рассчитанные значения состояний объекта и регулятора для времён из t. Состояния объекта содержатся в первых p.PLANT_DIM столбцах. В оставшихся c.HIDDEN_DIM столбцах содержатся состояния регулятора. |
u |
ndarray, shape=(t.shape[0], p.CONTROL_DIM) |
Двумерный массив, содержащий рассчитанные значения управления для времён из t. Возвращается, если return_control задано как True. |
animate_pendulum¶anim = animate_pendulum(t, y, xlim='auto', ylim=(-1.2, 1.2), resolution=(960, 540), dpi=108, spacing=1, invsec=1.0, filename=None, codec=None, progress=True)
Функция animate_pendulum строит объект-анимацию matplotlib.animation.Animation по заданным состояниям маятника на эквидистантной временной сетке. Длина маятника полагается равной 1.
Принимаемые аргументы
| Название | Тип | Описание | Значение по умолчанию |
|---|---|---|---|
t |
ndarray, ndim=1 |
Сортированный по возрастанию эквидистантный массив значений времени, для которых задаются состояния маятника. | |
y |
ndarray, ndim=2, shape[0]=t.shape[0], shape[1]>=3 |
Двумерный массив, содержащий состояния маятника для времён из t. По первому измерению меняется время, по второму — номер переменной, фактически используются только первый и третий столбцы, в которых содержатся значения угла $\varphi$ и положения каретки $\xi$ соответственно. |
|
xlim |
'auto' or 'center' or arraylike, shape=(2,) |
Интервал изменения горизонтальной координаты в анимации. Если xlim задано как 'auto' выбирается автоматически, в соответствии с заданным разрешением кадра и диапазоном движения каретки. Если xlim задано как 'center' выбирается автоматически, в соответствии с разрешением кадра, так что бы середина интервала совпдала с начальным положением каретки. |
'auto' |
ylim |
arraylike, shape=(2,) |
Интервал изменения вертикальной координаты в анимации. | (-1.2, 1.2) |
resolution |
arraylike, shape=(2,) |
Разрешение кадра. | (960, 540) |
dpi |
number |
Число пикселей на дюйм (влияет на толщину линий и размер символов). | 100 |
spacing |
int |
Количество шагов временной сетки t на один кадр анимации. |
1 |
invsec |
number |
Продолжительность единицы времени (в которой задана сетка t) в секундах. |
1.0 |
filename |
str or None |
Имя файла, в который записывается видео анимации, например 'test.mp4'. Если filename задано как None, запись в файл не производится. |
None |
codec |
str or None |
Кодировщик видео, например 'h264'. Если задано как None используется кодировщик по умолчанию. |
None |
progress |
bool |
Если задано как True, при сохранение в файл показывается индикатор прогресса. |
True |
Возращаемые значения
| Название | Тип | Описание |
|---|---|---|
anim |
matplotlib.animation.Animation |
Анимированный объект matplotlib. |
from control_theory import *
pd = Pendulum()
c0 = ZeroControl(4, 1) # Создание нулевого управления
td, yd = integrate(pd, c0, x0=[0., 0.1, 0., 0.], v0=empty(0), dt=0.01, T=30,
N=1)
plot(td, yd[:, 0])
xlabel(r'Время, $t$')
ylabel(r'Угол, $\varphi$')
title('Зависимость угла поворота от времени для маятника без управления');
subplot(2, 1, 1) # Первая панель (первый «подрисунок»)
plot(td, yd[:, 0])
xlabel(r'Время, $t$')
ylabel(r'Угол, $\varphi$')
subplot(2, 1, 2) # Вторая панель (второй «подрисунок»)
plot(td, yd[:, 2])
xlabel(r'Время, $t$')
ylabel(r'Положение каретки, $\xi$')
suptitle(r'Зависимости обобщённых координат от времени')
tight_layout(rect=[0, 0, 1, 0.95]); # Корректировка размеров осей с тем, чтобы
# подписи по осям не перекрывали других
# подрисунков и не выходили за пределы
# всего рисунка. Аргумент rect указывает
# на положение граничной рамки, содержащей
# все подрисунки, но не главный заголовок.
plot(yd[:, 2] + sin(yd[:, 0]), cos(yd[:, 0]))
gca().set_aspect('equal') # Установка одинаковых масштабов по вертикальной
# и горизонтальной осям
xlabel(r'$x$')
ylabel(r'$y$')
title('Траектории груза для маятника без управления');
# Создание анимации и запись в файл `test.mp4`
animd = animate_pendulum(td, yd, spacing=3, invsec=1,
filename='files/test.mp4')
# Импорт функции `HTML` для отображения произвольного кода HTML в выводе
# ячейки
from IPython.display import HTML
HTML(animd.to_html5_video()) # Вывод анимации в блокнот как встроенное видео
# HTML 5
rcParams['animation.embed_limit'] = 2.5 # Максимальный размер анимации в МБ
# Создание анимации меньщего размера
animd_small = animate_pendulum(td, yd, spacing=30, invsec=1)
# Вывод анимации в виде кода HTML
HTML(animd_small.to_jshtml())
В качестве альтернативы можно выводить видео из предварительно созданного файла в ячейках Markdown. В таком случае видео будет отображаться и после перезагрузки страницы. Видео ниже вставлено с помощью команды
<video controls frameborder="0" allowfullscreen autoplay src="test.mp4" type=video/mp4></video>integrate для расчёта эволюции демонстрационного маятника с управлением по углу¶c = LinearStateControl(4, 1, [2., 0., 0., 0.])
td, yd = integrate(pd, c, x0=[0., 0.1, 0., 0.], v0=empty(0), dt=0.01, T=30,
N=1)
plot(td, yd[:, 0])
xlabel(r'Время, $t$')
ylabel(r'Угол, $\varphi$')
title('Зависимость угла поворота от времени для маятника с управлением по ' +
'углу');
В системе $\dot x = y$, $\dot y = -x - y + x^3$ имеется две седловые точки $(x, y) = (\pm 1, 0)$. Система симметрична относительно замены $(x, y) \rightarrow (-x, -y)$, и можно рассмотреть только седло $(x, y) = (1, 0)$ и две его устойчивые сепаратрисы. Вблизи этого седла система линеаризуется как $\dot x = y$, $\dot y = (x - 1) - y$ или $$\frac{d}{dt} \begin{bmatrix} x - 1 \\ y \end{bmatrix} = \begin{bmatrix} 0 & 1 \\ 1 & -1 \end{bmatrix} \begin{bmatrix} x - 1 \\ y \end{bmatrix}.$$ Матрица системы имеет собственные значения $\dfrac{-1 \pm \sqrt{5}}{2}$ и собственные векторы $\begin{bmatrix} 1 \\ \frac{-1 \pm \sqrt{5}}{2} \end{bmatrix}$. Устойчивым сепаратрисам отвечает отрицательное собственное число (нижний знак). Таким образом, в качестве начальных условий для нахождения сепаратрис можно взять $$ \begin{bmatrix} x(0) - 1 \\ y(0) \end{bmatrix} = \pm\delta \begin{bmatrix} 1 \\ \frac{1 - \sqrt{5}}{2} \end{bmatrix}, $$ где $\delta$ — некоторое достаточное малое (в масштабах нелинейности системы, то есть единицы) положительное число.
delta = 1e-3 # Малое число
# Пределы по осям на рисунке
xmin = -3.
xmax = 3.
ymin = -3.
ymax = 3.
# Массив времён для сепаратрисы, приходящей в точку (1, 0) снизу справа
t1 = arange(0., -8.5, -0.01)
# Массив времён для сепаратрисы, приходящей в точку (1, 0) сверху слева
t2 = arange(0., -9.5, -0.01)
# Функция вычисляет правую часть системы (векторное поле потока)
def rhs_example(y, t):
return array([y[1], -y[1] - y[0] * (1. - y[0] * y[0])])
# Функция задаёт критерий остановки расчёта сепаратрис: расчёт
# останавливается, когда сепратриса покидает заданные границы осевой коробки
# (пределы по осям). В общем случае сепартриса, покинув заданные пределы,
# может затем вернуться обратно и могут потребоваться более сложные критерии
# остановки
def stop_condition_example(y, t):
return y[0] > xmax or y[0] < xmin or y[1] > ymax or y[1] < ymin
# Интегрирование методом Рунге — Кутты вдоль сепаратрис системы, определяемой
# функцией rhs_example в обратном времени. Для маятника можно вместо функции
# rk использовать функцию integrate, при этом нужно указать отрицательные
# значения шага по времени dt и времени расчёта T.
_, y1 = rk(rhs_example, [1. + delta, -0.5 * (1. + sqrt(5.)) * delta], t1, 10,
stop_condition_example)
_, y2 = rk(rhs_example, [1. - delta, 0.5 * (1. + sqrt(5.)) * delta], t2, 10,
stop_condition_example)
# Объединение двух сепаратрис в одну кривую
y = vstack([y1[::-1], y2])
# Объединение двух симметричных пар сепаратрис в одну кривую
y = vstack([y, -y])
# Заполнение области, ограниченной сепаратрисами цветом
fill(y[:, 0], y[:, 1], color='pink')
# Построение сепаратрис
plot(y[:, 0], y[:, 1], color='red')
# Точки в положениях равновесия
plot([0], [0], 'bo')
plot([-1, 1], [0, 0], 'ro')
# Установление пределов по осям
axis([xmin, xmax, ymin, ymax])
# Построение фазовой плоскости
x, y = meshgrid(linspace(xmin, xmax, 201), linspace(ymin, ymax, 201))
dy = -y - x * (1. - x * x)
streamplot(x, y, y, dy, color='k', density=1.5)
xlabel(r'$x$')
ylabel(r'$y$')
title(r'Область притяжения устойчивого положения равновесия системы ' +
r'$\dot x = y$, $\dot y = -x - y + x^3$');
Характеристический полином уравнения $(9)$ имеет вид $Q(p) = (1 - \alpha) p^2 + (\lambda + k_1) p - \Omega^2 + k_0$. Параметризуя границу интересующей области как $p = i\omega$, получаем систему уравнений $$k_0 - \Omega^2 - (1 - \alpha) \omega^2 = 0,$$ $$\omega k_1 + \omega \lambda = 0.$$ Уравнение особой прямой при $\omega = 0$ имеет вид $$k_0 - \Omega^2 = 0.$$ Решение системы при $\omega \neq 0$ имеет вид $$k_0 = \omega^2 (1 - \alpha) + \Omega^2,$$ $$k_1 = -\lambda$$ и определяет проходимый дважды луч с началом в $(k_0, k_1) = (-\Omega^2, -\lambda)$ и направлением вдоль положительного направления оси $k_0$.
# Инициализация маятника без трения о внешнюю среду
pd2 = Pendulum()
pd2.mu = 0.
pd2.nu = 0.
k0lim = (-5., 5.) # Пределы по горизонтальной оси
k1lim = (-5., 5.) # Пределы по вертикальной оси
# Две точки на особой прямой
k1_0 = array(k1lim) # ᅟᅟОрдинаты
k0_0 = full((2,), pd2.omega * pd2.omega) # Абсциссы
# Две точки на горизонтальном луче
k0_1 = array([pd2.omega * pd2.omega, k0lim[1]]) # Абсциссы
k1_1 = full((2,), -pd2.lambd) # ᅟОрдинаты
plot(k0_0, k1_0, k0_1, k1_1)
xlim(k0lim)
ylim(k1lim)
xlabel(r'$k_0$')
ylabel(r'$k_1$')
text(-2.5, 0., '1')
text(2.5, 2.5, '2')
text(2.5, -2.5, '0')
title(r'D-разбиение для коэффициентов управления по углу и угловой ' +
r'скорости, $\mathrm{Re}\, p < 0$');
Pc = 0.1 * arange(-10, 1) # Значения p_c, для которых будут построены
# линии уровня
alim = (-1., 5.) # Пределы по горизонтальной оси
blim = (-1., 5.) # Пределы по вертикальной оси
na = 201 # Количество шагов сетки по горизонтали
nb = 201 # Количество шагов сетки по вертикали
# Создание сеток (матриц) значений a и b
a = convolve(linspace(alim[0], alim[1], na), [0.5, 0.5], 'valid')
b = convolve(linspace(blim[0], blim[1], nb), [0.5, 0.5], 'valid')
A, B = meshgrid(a, b)
# Функция находит максимальную действительную часть корней полинома по
# заданным значениям a и b
def find_maximal_re(a0, b0):
return amax(real(roots([1., a0, b0])))
# Векторизация функции find_maximal_re для работы с массивами значений a и b
vfind_maximal_re = vectorize(find_maximal_re)
# Вычисление матрицы значений максимальной действительной части нуля полинома
maxrez = vfind_maximal_re(A, B)
Контурную карту с метками линий уровня можно построить с помощью функций contour и clabel.
contours = contour(a, b, maxrez, Pc)
clabel(contours, inline=1, fontsize=8)
xlabel(r'$a$')
ylabel(r'$b$')
title(r'Границы областей $\max \,\mathrm{Re}\, p < p_c$ для уравнения ' +
r'$p^2 + ap + b = 0$');
Вместо меток линий уровня можно отобразить шкалу или построить соответствующую тепловую карту со шкалой с помощью функций imshow и colorbar.
im = imshow(maxrez, interpolation='bilinear', origin='lower', cmap=cm.gray,
aspect='auto', extent=(alim[0], alim[1], blim[0], blim[1]))
contours = contour(maxrez, Pc, origin='lower', linewidths=2,
extent=(alim[0], alim[1], blim[0], blim[1]))
xlabel(r'$a$')
ylabel(r'$b$')
title(r'Тепловая карта $\max \,\mathrm{Re}\, p$ для уравнения ' +
r'$p^2 + ap + b = 0$')
contours_bar = colorbar(contours, orientation='vertical', shrink=0.8,
extend='both', label=r'$\max \,\mathrm{Re}\, p$')
im_bar = colorbar(im, orientation='vertical', shrink=0.8,
label=r'$\max \,\mathrm{Re}\, p$');
p = Pendulum()
c = LinearSinPendulumControl()
t, y = integrate(p, c, x0=[0., 0.1, 0., 0.], v0=empty(0), dt=0.01, T=30, N=1)
animate_pendulum(t, y, resolution=(3840, 2160), dpi=432, spacing=3, invsec=1,
filename='files/example2160p.mp4')
animate_pendulum(t, y, resolution=(1920, 1080), dpi=216, spacing=3, invsec=1,
filename='files/example1080p.mp4')
animate_pendulum(t, y, resolution=(1280, 720), dpi=144, spacing=3, invsec=1,
filename='files/example720p.mp4')
animate_pendulum(t, y, resolution=(960, 540), dpi=108, spacing=3, invsec=1,
filename='files/example540p.mp4')
animate_pendulum(t, y, resolution=(854, 480), dpi=96, spacing=3, invsec=1,
filename='files/example480p.mp4')
animate_pendulum(t, y, resolution=(640, 360), dpi=72, spacing=3, invsec=1,
filename='files/example360p.mp4');
# Функции для определения бифуркационных значений c.k[1]
c = LinearSinPendulumControl([2, 0.3685484, 0, 0])
delta = 1e-4 # Малое число, характеризующее близость стартовой точки
# интегрирования к седловой точке
dt = -0.01 # Шаг интегрирования
tol = 0.001 # Погрешность в определении gamma_1 и gamma_2
T = -40 # Конечное время для интегрирования
powi = 2
imax = 30
def integrate_for_gamma1(k1):
c.k[1] = k1
tphi = arccos(pd2.omega * pd2.omega / c.k[0])
pminus = -(c.k[0] * (pd2.lambd * c.k[0] + pd2.omega ** 2 * c.k[1]) +
sqrt((pd2.lambd * c.k[0] + pd2.omega ** 2 * c.k[1]) ** 2 *
c.k[0] ** 2 + 4. * c.k[0] * (c.k[0] ** 2 - pd2.alpha *
pd2.omega ** 4) *
(c.k[0] ** 2 - pd2.omega ** 4))) / (2 * (c.k[0] ** 2 -
pd2.alpha *
pd2.omega ** 4))
x0 = [tphi + delta, delta * pminus, 0., 0.]
_, y = integrate(pd2, c, x0=x0, v0=empty(0), dt=dt, T=T, N=1,
stop_condition=lambda y, t: y[0] > 2 * pi - tphi or
y[1] > 0)
return y[-1]
def find_gamma1(k0, k1est):
c.k[0] = k0
tphi = arccos(pd2.omega * pd2.omega / c.k[0])
z0 = integrate_for_gamma1(k1est)
k1start = k1est
k1end = k1est
if z0[0] < 2 * pi - tphi:
i = 1
while z0[0] < 2 * pi - tphi:
k1start = k1end - (powi ** i) * tol
z0 = integrate_for_gamma1(k1start)
i += 1
if i == imax:
print('Error: find_gamma1, imax too low!')
print('k0 = {0}, k1est = {1}'.format(k0, k1est))
return None
else:
i = 1
while z0[0] >= 2 * pi - tphi:
k1end = k1start + (powi ** i) * tol
z0 = integrate_for_gamma1(k1end)
i += 1
if i == imax:
print('Error: find_gamma1, imax too low!')
print('k0 = {0}, k1est = {1}'.format(k0, k1est))
return None
while k1end - k1start > tol:
z0 = integrate_for_gamma1(0.5 * (k1start + k1end))
if z0[0] < 2 * pi - tphi:
k1end = 0.5 * (k1start + k1end)
else:
k1start = 0.5 * (k1start + k1end)
return 0.5 * (k1start + k1end)
def integrate_for_gamma2(k1):
c.k[1] = k1
tphi = arccos(pd2.omega * pd2.omega / c.k[0])
pminus = -(c.k[0] * (pd2.lambd * c.k[0] + pd2.omega ** 2 * c.k[1]) +
sqrt((pd2.lambd * c.k[0] + pd2.omega ** 2 * c.k[1]) ** 2 *
c.k[0] ** 2 + 4. * c.k[0] * (c.k[0] ** 2 - pd2.alpha *
pd2.omega ** 4) *
(c.k[0] ** 2 - pd2.omega ** 4))) / (2 * (c.k[0] ** 2 -
pd2.alpha *
pd2.omega ** 4))
x0 = [tphi - delta, -delta * pminus, 0., 0.]
_, y = integrate(pd2, c, x0=x0, v0=empty(0), dt=dt, T=T, N=1,
stop_condition=lambda y, t: y[0] < tphi - 2 * pi or
y[1] < 0)
return y[-1]
def find_gamma2(k0, k1est):
c.k[0] = k0
tphi = arccos(pd2.omega * pd2.omega / c.k[0])
z0 = integrate_for_gamma2(k1est)
k1start = k1est
k1end = k1est
if z0[0] > tphi - 2 * pi:
i = 1
while z0[0] > tphi - 2 * pi:
k1start = k1end - (powi ** i) * tol
z0 = integrate_for_gamma2(k1start)
i += 1
if i == imax:
print('Error: find_gamma2, imax too low!')
print('k0 = {0}, k1est = {1}'.format(k0, k1est))
return None
else:
i = 1
while z0[0] <= tphi - 2 * pi:
k1end = k1start + (powi ** i) * tol
z0 = integrate_for_gamma2(k1end)
i += 1
if i == imax:
print('Error: find_gamma2, imax too low!')
print('k0 = {0}, k1est = {1}'.format(k0, k1est))
return None
while k1end - k1start > tol:
z0 = integrate_for_gamma2(0.5 * (k1start + k1end))
if z0[0] > tphi - 2 * pi:
k1end = 0.5 * (k1start + k1end)
else:
k1start = 0.5 * (k1start + k1end)
return 0.5 * (k1start + k1end)
# Бифуркационные значения для c.k[0] = 2
tol = 1e-7
print('gamma1 = {0:.8g}; gamma2 = {1:.8g}'.format(find_gamma1(2., 0.2),
find_gamma2(2., 0.3)))
print('lambda = {0:.8g}'.format(-pd2.lambd))
gamma1 = 0.19656405; gamma2 = 0.46080815 lambda = -0.13476466
c = LinearSinPendulumControl([2., 0., 0., 0.])
tphi = arccos(pd2.omega * pd2.omega / c.k[0])
pminus = -(c.k[0] * (pd2.lambd * c.k[0] + pd2.omega ** 2 * c.k[1]) +
sqrt((pd2.lambd * c.k[0] + pd2.omega ** 2 * c.k[1]) ** 2 *
c.k[0] ** 2 + 4. * c.k[0] * (c.k[0] ** 2 -
pd2.alpha * pd2.omega ** 4) *
(c.k[0] ** 2 - pd2.omega ** 4))) / (2 * (c.k[0] ** 2 -
pd2.alpha *
pd2.omega ** 4))
delta = 1e-4
phimin = -pi
phimax = pi
dphimin = -10.
dphimax = 10.
dt1 = -0.01
dt2 = -0.01
T1 = -20.
T2 = -20.
# Критерии остановки
def stop_condition1(y, t):
if y[1] >= dphimin:
stop_condition1.philast = y[0]
return False
elif y[0] < stop_condition1.philast + 2 * pi:
return False
else:
return True
def stop_condition2(y, t):
if y[1] <= dphimax:
stop_condition2.philast = y[0]
return False
elif y[0] > stop_condition2.philast - 2 * pi:
return False
else:
return True
# ᅟРасчёт сепаратрис
stop_condition1.philast = tphi
_, y1 = integrate(pd2, c, x0=[tphi + delta, delta * pminus, 0., 0.],
v0=empty(0), dt=dt1, T=T1, N=10,
stop_condition=stop_condition1)
stop_condition2.philast = tphi
_, y2 = integrate(pd2, c, x0=[tphi - delta, -delta * pminus, 0., 0.],
v0=empty(0), dt=dt2, T=T2, N=10,
stop_condition=stop_condition2)
# Объединение двух сепаратрис в одну кривую
y = vstack([y1[::-1], y2])
# Объединение двух симметричных пар сепаратрис в одну кривую
y = vstack([y, -y])
# Заполнение области, ограниченной сепаратрисами цветом
i1 = int(floor((amin(y[:, 0]) + pi) / (2 * pi)))
i2 = int(ceil((amax(y[:, 0]) - pi) / (2 * pi)))
for i in range(i1, i2 + 1):
fill(y[:, 0] - 2 * i * pi, y[:, 1], color='pink')
# Построение сепаратрис
for i in range(i1, i2 + 1):
plot(y[:, 0] - 2 * i * pi, y[:, 1], color='red')
# Точки в положениях равновесия
plot([0, -pi, pi], [0, 0, 0], 'bo')
plot([-tphi, tphi], [0, 0], 'ro')
# Установление пределов по осям
axis([phimin, phimax, dphimin, dphimax])
# Построение фазовой плоскости
phi, dphi = meshgrid(linspace(phimin, phimax, 201),
linspace(dphimin, dphimax, 201))
ddphi = ((pd2.omega * pd2.omega - (c.k[0] + pd2.alpha * dphi * dphi) *
cos(phi)) * sin(phi) - (pd2.lambd + c.k[1] * cos(phi)) *
dphi) / (1. - pd2.alpha * cos(phi) * cos(phi))
streamplot(phi, dphi, dphi, ddphi, color='k', density=2., linewidth=0.5)
xlabel(r'$\varphi$')
ylabel(r'$\dot\varphi$')
title(r'Фазовый портрет замкнутой системы при $k_0 > \Omega^2$, ' +
r'$-\lambda < k_1 < \lambda$');
savefig('files/phaseplane1.png', dpi='figure', format='png', transparent=True,
frameon=False)
savefig('files/phaseplane1.pdf', dpi='figure', format='pdf', transparent=True,
frameon=False)
c = LinearSinPendulumControl([2, 0.15, 0, 0])
tphi = arccos(pd2.omega * pd2.omega / c.k[0])
pminus = -(c.k[0] * (pd2.lambd * c.k[0] + pd2.omega ** 2 * c.k[1]) +
sqrt((pd2.lambd * c.k[0] + pd2.omega ** 2 * c.k[1]) ** 2 *
c.k[0] ** 2 + 4. * c.k[0] * (c.k[0] ** 2 -
pd2.alpha * pd2.omega ** 4) *
(c.k[0] ** 2 - pd2.omega ** 4))) / (2 * (c.k[0] ** 2 -
pd2.alpha *
pd2.omega ** 4))
delta = 1e-4
phimin = -pi
phimax = pi
dphimin = -10.
dphimax = 10.
dt1 = -0.01
dt2 = -0.01
T1 = -20.
T2 = -20.
# Критерии остановки
def stop_condition1(y, t):
if y[1] >= dphimin:
stop_condition1.philast = y[0]
return False
elif y[0] < stop_condition1.philast + 2 * pi:
return False
else:
return True
def stop_condition2(y, t):
if y[1] <= dphimax:
stop_condition2.philast = y[0]
return False
elif y[0] > stop_condition2.philast - 2 * pi:
return False
else:
return True
# ᅟРасчёт сепаратрис
stop_condition1.philast = tphi
_, y1 = integrate(pd2, c, x0=[tphi + delta, delta * pminus, 0., 0.],
v0=empty(0), dt=dt1, T=T1, N=10,
stop_condition=stop_condition1)
stop_condition2.philast = tphi
_, y2 = integrate(pd2, c, x0=[tphi - delta, -delta * pminus, 0., 0.],
v0=empty(0), dt=dt2, T=T2, N=10,
stop_condition=stop_condition2)
# Расчёт предельного цикла
_, y3 = integrate(pd2, c, x0=[4.07374665, 0.57169832, 0, 0], v0=empty(0),
dt=0.01, T=3.215, N=10)
y3 = vstack([y3, y3[-1]])
# Объединение двух сепаратрис в одну кривую
y = vstack([y1[::-1], y2])
# Объединение двух симметричных пар сепаратрис в одну кривую
y = vstack([y, -y])
# Заполнение области, ограниченной сепаратрисами цветом
i1 = int(floor((amin(y[:, 0]) + pi) / (2 * pi)))
i2 = int(ceil((amax(y[:, 0]) - pi) / (2 * pi)))
for i in range(i1, i2 + 1):
fill(y[:, 0] - 2 * i * pi, y[:, 1], color='pink')
# Построение сепаратрис
for i in range(i1, i2 + 1):
plot(y[:, 0] - 2 * i * pi, y[:, 1], color='red')
# Построение предельного цикла
i1 = int(floor((amin(y3[:, 0]) + pi) / (2 * pi)))
i2 = int(ceil((amax(y3[:, 0]) - pi) / (2 * pi)))
for i in range(i1, i2 + 1):
plot(y3[:, 0] - 2 * i * pi, y3[:, 1], color='blue')
# Точки в положениях равновесия
plot([0], [0], 'bo')
plot([-tphi, tphi, -pi, pi], [0, 0, 0, 0], 'ro')
# Установление пределов по осям
axis([phimin, phimax, dphimin, dphimax])
# Построение фазовой плоскости
phi, dphi = meshgrid(linspace(phimin, phimax, 201),
linspace(dphimin, dphimax, 201))
ddphi = ((pd2.omega * pd2.omega - (c.k[0] + pd2.alpha * dphi * dphi) *
cos(phi)) * sin(phi) - (pd2.lambd + c.k[1] * cos(phi)) *
dphi) / (1. - pd2.alpha * cos(phi) * cos(phi))
streamplot(phi, dphi, dphi, ddphi, color='k', density=2., linewidth=0.5)
xlabel(r'$\varphi$')
ylabel(r'$\dot\varphi$')
title(r'Фазовый портрет замкнутой системы при $k_0 > \Omega^2$, ' +
r'$\lambda < k_1 < \gamma_1(k_0)$');
savefig('files/phaseplane2.png', dpi='figure', format='png', transparent=True,
frameon=False)
savefig('files/phaseplane2.pdf', dpi='figure', format='pdf', transparent=True,
frameon=False)
c = LinearSinPendulumControl([2, 0.19656405, 0, 0])
tphi = arccos(pd2.omega * pd2.omega / c.k[0])
pminus = -(c.k[0] * (pd2.lambd * c.k[0] + pd2.omega ** 2 * c.k[1]) +
sqrt((pd2.lambd * c.k[0] + pd2.omega ** 2 * c.k[1]) ** 2 *
c.k[0] ** 2 + 4. * c.k[0] * (c.k[0] ** 2 -
pd2.alpha * pd2.omega ** 4) *
(c.k[0] ** 2 - pd2.omega ** 4))) / (2 * (c.k[0] ** 2 -
pd2.alpha *
pd2.omega ** 4))
delta = 1e-4
phimin = -pi
phimax = pi
dphimin = -10.
dphimax = 10.
dt1 = -0.01
dt2 = -0.01
T1 = -20.
T2 = -20.
# Критерии остановки
def stop_condition1(y, t):
return y[0] > 2 * pi - tphi or y[1] > 0
def stop_condition2(y, t):
if y[1] <= dphimax:
stop_condition2.philast = y[0]
return False
elif y[0] > stop_condition2.philast - 2 * pi:
return False
else:
return True
# Расчёт гетероклинической траектории
_, y1 = integrate(pd2, c, x0=[tphi + delta, delta * pminus, 0., 0.],
v0=empty(0), dt=dt1, T=T1, N=10,
stop_condition=stop_condition1)
# Расчёт сепаратрисы
stop_condition2.philast = tphi
_, y2 = integrate(pd2, c, x0=[tphi - delta, -delta * pminus, 0., 0.],
v0=empty(0), dt=dt2, T=T2, N=10,
stop_condition=stop_condition2)
# Объединение двух гетероклинических траекторий в один цикл
yphi1 = hstack([[tphi], y1[:, 0], [3 * pi - tphi], 2 * pi - y1[:, 0]])
ydphi1 = hstack([[0], y1[:, 1], [0], -y1[:, 1]])
# Заполнение всей плоскости розовым цветом
fill([phimin, phimax, phimax, phimin], [dphimin, dphimin, dphimax, dphimax],
color='pink')
# Заполнение внутренней области гетероклинического контура белым цветом
fill(yphi1, ydphi1, 'w', -yphi1, -ydphi1, 'w')
# Построение гетероклинического контура
plot(yphi1, ydphi1, 'r', -yphi1, -ydphi1, 'r')
# Построение оставшихся сепаратрис
i1 = int(floor((amin(y2[:, 0]) + pi) / (2 * pi)))
i2 = int(ceil((amax(-y2[:, 0]) - pi) / (2 * pi)))
for i in range(i1, i2 + 1):
plot(y2[:, 0] - 2 * i * pi, y2[:, 1], color='red')
plot(-y2[:, 0] - 2 * i * pi, -y2[:, 1], color='red')
# Точки в положениях равновесия
plot([0], [0], 'bo')
plot([-tphi, tphi, -pi, pi], [0, 0, 0, 0], 'ro')
# Установление пределов по осям
axis([phimin, phimax, dphimin, dphimax])
# Построение фазовой плоскости
phi, dphi = meshgrid(linspace(phimin, phimax, 201),
linspace(dphimin, dphimax, 201))
ddphi = ((pd2.omega * pd2.omega - (c.k[0] + pd2.alpha * dphi * dphi) *
cos(phi)) * sin(phi) - (pd2.lambd + c.k[1] * cos(phi)) *
dphi) / (1. - pd2.alpha * cos(phi) * cos(phi))
streamplot(phi, dphi, dphi, ddphi, color='k', density=2., linewidth=0.5)
xlabel(r'$\varphi$')
ylabel(r'$\dot\varphi$')
title(r'Фазовый портрет замкнутой системы при $k_0 > \Omega^2$, ' +
r'$k_1 = \gamma_1(k_0)$');
savefig('files/phaseplane3.png', dpi='figure', format='png', transparent=True,
frameon=False)
savefig('files/phaseplane3.pdf', dpi='figure', format='pdf', transparent=True,
frameon=False)
c = LinearSinPendulumControl([2, 0.3, 0, 0])
tphi = arccos(pd2.omega * pd2.omega / c.k[0])
pminus = -(c.k[0] * (pd2.lambd * c.k[0] + pd2.omega ** 2 * c.k[1]) +
sqrt((pd2.lambd * c.k[0] + pd2.omega ** 2 * c.k[1]) ** 2 *
c.k[0] ** 2 + 4. * c.k[0] * (c.k[0] ** 2 -
pd2.alpha * pd2.omega ** 4) *
(c.k[0] ** 2 - pd2.omega ** 4))) / (2 * (c.k[0] ** 2 -
pd2.alpha *
pd2.omega ** 4))
delta = 1e-4
phimin = -pi
phimax = pi
dphimin = -10.
dphimax = 10.
dt1 = -0.01
dt2 = -0.01
T1 = -50.
T2 = -20.
# Критерии остановки
def stop_condition1(y, t):
return (y[0] - pi) * (y[0] - pi) + y[1] * y[1] < delta * delta
def stop_condition2(y, t):
if y[1] <= dphimax:
stop_condition2.philast = y[0]
return False
elif y[0] > stop_condition2.philast - 2 * pi:
return False
else:
return True
# Расчёт сепаратрис
_, y1 = integrate(pd2, c, x0=[tphi + delta, delta * pminus, 0., 0.],
v0=empty(0), dt=dt1, T=T1, N=10,
stop_condition=stop_condition1)
stop_condition2.philast = tphi
_, y2 = integrate(pd2, c, x0=[tphi - delta, -delta * pminus, 0., 0.],
v0=empty(0), dt=dt2, T=T2, N=10,
stop_condition=stop_condition2)
# Объединение двух сепаратрис в одну кривую
y = vstack([y1[::-1], y2])
# Заполнение всей плоскости розовым цветом
fill([phimin, phimax, phimax, phimin], [dphimin, dphimin, dphimax, dphimax],
color='pink')
# Построение сепаратрис
i1 = int(floor((amin([y[:, 0], -y[:, 0]]) + pi) / (2 * pi)))
i2 = int(ceil((amax([y[:, 0], -y[:, 0]]) - pi) / (2 * pi)))
for i in range(i1, i2 + 1):
plot(y[:, 0] - 2 * i * pi, y[:, 1], color='red')
plot(-y[:, 0] - 2 * i * pi, -y[:, 1], color='red')
# Точки в положениях равновесия
plot([0], [0], 'bo')
plot([-tphi, tphi, -pi, pi], [0, 0, 0, 0], 'ro')
# Установление пределов по осям
axis([phimin, phimax, dphimin, dphimax])
# Построение фазовой плоскости
phi, dphi = meshgrid(linspace(phimin, phimax, 201),
linspace(dphimin, dphimax, 201))
ddphi = ((pd2.omega * pd2.omega - (c.k[0] + pd2.alpha * dphi * dphi) *
cos(phi)) * sin(phi) - (pd2.lambd + c.k[1] * cos(phi)) *
dphi) / (1. - pd2.alpha * cos(phi) * cos(phi))
streamplot(phi, dphi, dphi, ddphi, color='k', density=2., linewidth=0.5)
xlabel(r'$\varphi$')
ylabel(r'$\dot\varphi$')
title(r'Фазовый портрет замкнутой системы при $k_0 > \Omega^2$, ' +
r'$\gamma_1(k_0) < k_1 < \gamma_2(k_0)$');
savefig('files/phaseplane4.png', dpi='figure', format='png', transparent=True,
frameon=False)
savefig('files/phaseplane4.pdf', dpi='figure', format='pdf', transparent=True,
frameon=False)
c = LinearSinPendulumControl([2, 0.46080815, 0, 0])
tphi = arccos(pd2.omega * pd2.omega / c.k[0])
pminus = -(c.k[0] * (pd2.lambd * c.k[0] + pd2.omega ** 2 * c.k[1]) +
sqrt((pd2.lambd * c.k[0] + pd2.omega ** 2 * c.k[1]) ** 2 *
c.k[0] ** 2 + 4. * c.k[0] * (c.k[0] ** 2 -
pd2.alpha * pd2.omega ** 4) *
(c.k[0] ** 2 - pd2.omega ** 4))) / (2 * (c.k[0] ** 2 -
pd2.alpha *
pd2.omega ** 4))
delta = 1e-4
phimin = -pi
phimax = pi
dphimin = -10.
dphimax = 10.
dt1 = -0.01
dt2 = -0.01
T1 = -50.
T2 = -20.
# Критерии остановки
def stop_condition1(y, t):
return (y[0] - pi) * (y[0] - pi) + y[1] * y[1] < delta * delta
def stop_condition2(y, t):
return y[0] < tphi - 2 * pi or y[1] < 0
# Расчёт сепаратрисы
_, y1 = integrate(pd2, c, x0=[tphi + delta, delta * pminus, 0., 0.],
v0=empty(0), dt=dt1, T=T1, N=10,
stop_condition=stop_condition1)
# Расчёт гомоклинической траектории
_, y2 = integrate(pd2, c, x0=[tphi - delta, -delta * pminus, 0., 0.],
v0=empty(0), dt=dt2, T=T2, N=10,
stop_condition=stop_condition2)
# Объединение двух гомоклинических траекторий
yphi2 = hstack([y2[:, 0] + 2 * pi, [tphi], y2[:, 0]])
ydphi2 = hstack([y2[:, 1], [0], y2[:, 1]])
yphi2pm = hstack([yphi2, -yphi2])
ydphi2pm = hstack([ydphi2, -ydphi2])
# Заполнение цветом области между гомоклиническими траекториями
fill(yphi2pm, ydphi2pm, color='pink')
# Построение гомоклинических траекторий
plot(yphi2pm, ydphi2pm, color='red')
# Построение оставшихся сепаратрис
i1 = int(floor((amin([y1[:, 0], -y1[:, 0]]) + pi) / (2 * pi)))
i2 = int(ceil((amax([y1[:, 0], -y1[:, 0]]) - pi) / (2 * pi)))
for i in range(i1, i2 + 1):
plot(y1[:, 0] - 2 * i * pi, y1[:, 1], color='red')
plot(-y1[:, 0] - 2 * i * pi, -y1[:, 1], color='red')
# Точки в положениях равновесия
plot([0], [0], 'bo')
plot([-tphi, tphi, -pi, pi], [0, 0, 0, 0], 'ro')
# Установление пределов по осям
axis([phimin, phimax, dphimin, dphimax])
# Построение фазовой плоскости
phi, dphi = meshgrid(linspace(phimin, phimax, 201),
linspace(dphimin, dphimax, 201))
ddphi = ((pd2.omega * pd2.omega - (c.k[0] + pd2.alpha * dphi * dphi) *
cos(phi)) * sin(phi) - (pd2.lambd + c.k[1] * cos(phi)) *
dphi) / (1. - pd2.alpha * cos(phi) * cos(phi))
streamplot(phi, dphi, dphi, ddphi, color='k', density=2., linewidth=0.5)
xlabel(r'$\varphi$')
ylabel(r'$\dot\varphi$')
title(r'Фазовый портрет замкнутой системы при $k_0 > \Omega^2$, ' +
r'$k_1 = \gamma_2(k_0)$');
savefig('files/phaseplane5.png', dpi='figure', format='png', transparent=True,
frameon=False)
savefig('files/phaseplane5.pdf', dpi='figure', format='pdf', transparent=True,
frameon=False)
c = LinearSinPendulumControl([2, 0.6, 0, 0])
tphi = arccos(pd2.omega * pd2.omega / c.k[0])
pminus = -(c.k[0] * (pd2.lambd * c.k[0] + pd2.omega ** 2 * c.k[1]) +
sqrt((pd2.lambd * c.k[0] + pd2.omega ** 2 * c.k[1]) ** 2 *
c.k[0] ** 2 + 4. * c.k[0] * (c.k[0] ** 2 -
pd2.alpha * pd2.omega ** 4) *
(c.k[0] ** 2 - pd2.omega ** 4))) / (2 * (c.k[0] ** 2 -
pd2.alpha *
pd2.omega ** 4))
delta = 1e-4
phimin = -pi
phimax = pi
dphimin = -10.
dphimax = 10.
dt1 = -0.01
dt2 = -0.01
T1 = -50.
T2 = -50.
# Критерии остановки
def stop_condition1(y, t):
return (y[0] - pi) * (y[0] - pi) + y[1] * y[1] < delta * delta
def stop_condition2(y, t):
return (y[0] + pi) * (y[0] + pi) + y[1] * y[1] < delta * delta
# Расчёт сепаратрис
_, y1 = integrate(pd2, c, x0=[tphi + delta, delta * pminus, 0., 0.],
v0=empty(0), dt=dt1, T=T1, N=10,
stop_condition=stop_condition1)
_, y2 = integrate(pd2, c, x0=[tphi - delta, -delta * pminus, 0., 0.],
v0=empty(0), dt=dt2, T=T2, N=10,
stop_condition=stop_condition2)
# Расчёт предельного цикла
_, y3 = integrate(pd2, c, x0=[-0.0760595019, 2.5430062513, 0, 0], v0=empty(0),
dt=0.01, T=4.205, N=10)
# Объединение двух сепаратрис в одну кривую
y = vstack([y1[::-1], y2])
# Объединение двух симметричных пар сепаратрис в одну кривую
y = vstack([y, -y])
# Объединение двух периодов предельного цикла в одну кривую
yphi3 = hstack([y3[:, 0] - 2 * pi, y3[:, 0]])
ydphi3 = hstack([y3[:, 1], y3[:, 1]])
# Заполнение области, ограниченной сепаратрисами цветом
i1 = int(floor((amin(y[:, 0]) + pi) / (2 * pi)))
i2 = int(ceil((amax(y[:, 0]) - pi) / (2 * pi)))
for i in range(i1, i2 + 1):
fill(y[:, 0] - 2 * i * pi, y[:, 1], color='pink')
# Построение сепаратрис
for i in range(i1, i2 + 1):
plot(y[:, 0] - 2 * i * pi, y[:, 1], color='red')
# Построение предельного цикла
plot(yphi3, ydphi3, color='blue')
plot(-yphi3, -ydphi3, color='blue')
# Точки в положениях равновесия
plot([0], [0], 'bo')
plot([-tphi, tphi, -pi, pi], [0, 0, 0, 0], 'ro')
# Установление пределов по осям
axis([phimin, phimax, dphimin, dphimax])
# Построение фазовой плоскости
phi, dphi = meshgrid(linspace(phimin, phimax, 201),
linspace(dphimin, dphimax, 201))
ddphi = ((pd2.omega * pd2.omega - (c.k[0] + pd2.alpha * dphi * dphi) *
cos(phi)) * sin(phi) - (pd2.lambd + c.k[1] * cos(phi)) *
dphi) / (1. - pd2.alpha * cos(phi) * cos(phi))
streamplot(phi, dphi, dphi, ddphi, color='k', density=2., linewidth=0.5)
xlabel(r'$\varphi$')
ylabel(r'$\dot\varphi$')
title(r'Фазовый портрет замкнутой системы при $k_0 > \Omega^2$, ' +
r'$k_1 > \gamma_2(k_0)$');
savefig('files/phaseplane6.png', dpi='figure', format='png', transparent=True,
frameon=False)
savefig('files/phaseplane6.pdf', dpi='figure', format='pdf', transparent=True,
frameon=False)
# Функции для построения бифуркационных диаграмм
c = LinearSinPendulumControl([2, 0.3, 0, 0])
delta = 1e-4 # Малое число, характеризующее близость стартовой точки
# интегрирования к седловой точке
dt = -0.01 # Шаг интегрирования
tol = 0.001 # Погрешность в определении gamma_1 и gamma_2
T = -40 # Конечное время для интегрирования
powi = 2
imax = 30
def integrate_for_gamma1(k1):
c.k[1] = k1
tphi = arccos(pd2.omega * pd2.omega / c.k[0])
pminus = -(c.k[0] * (pd2.lambd * c.k[0] + pd2.omega ** 2 * c.k[1]) +
sqrt((pd2.lambd * c.k[0] + pd2.omega ** 2 * c.k[1]) ** 2 *
c.k[0] ** 2 + 4. * c.k[0] * (c.k[0] ** 2 - pd2.alpha *
pd2.omega ** 4) *
(c.k[0] ** 2 - pd2.omega ** 4))) / (2 * (c.k[0] ** 2 -
pd2.alpha *
pd2.omega ** 4))
x0 = [tphi + delta, delta * pminus, 0., 0.]
_, y = integrate(pd2, c, x0=x0, v0=empty(0), dt=dt, T=T, N=1,
stop_condition=lambda y, t: y[0] > 2 * pi - tphi or
y[1] > 0)
return y[-1]
def find_gamma1(k0, k1est):
c.k[0] = k0
tphi = arccos(pd2.omega * pd2.omega / c.k[0])
z0 = integrate_for_gamma1(k1est)
k1start = k1est
k1end = k1est
if z0[0] < 2 * pi - tphi:
i = 1
while z0[0] < 2 * pi - tphi:
k1start = k1end - (powi ** i) * tol
z0 = integrate_for_gamma1(k1start)
i += 1
if i == imax:
print('Error: find_gamma1, imax too low!')
print('k0 = {0}, k1est = {1}'.format(k0, k1est))
return None
else:
i = 1
while z0[0] >= 2 * pi - tphi:
k1end = k1start + (powi ** i) * tol
z0 = integrate_for_gamma1(k1end)
i += 1
if i == imax:
print('Error: find_gamma1, imax too low!')
print('k0 = {0}, k1est = {1}'.format(k0, k1est))
return None
while k1end - k1start > tol:
z0 = integrate_for_gamma1(0.5 * (k1start + k1end))
if z0[0] < 2 * pi - tphi:
k1end = 0.5 * (k1start + k1end)
else:
k1start = 0.5 * (k1start + k1end)
return 0.5 * (k1start + k1end)
def integrate_for_gamma2(k1):
c.k[1] = k1
tphi = arccos(pd2.omega * pd2.omega / c.k[0])
pminus = -(c.k[0] * (pd2.lambd * c.k[0] + pd2.omega ** 2 * c.k[1]) +
sqrt((pd2.lambd * c.k[0] + pd2.omega ** 2 * c.k[1]) ** 2 *
c.k[0] ** 2 + 4. * c.k[0] * (c.k[0] ** 2 - pd2.alpha *
pd2.omega ** 4) *
(c.k[0] ** 2 - pd2.omega ** 4))) / (2 * (c.k[0] ** 2 -
pd2.alpha *
pd2.omega ** 4))
x0 = [tphi - delta, -delta * pminus, 0., 0.]
_, y = integrate(pd2, c, x0=x0, v0=empty(0), dt=dt, T=T, N=1,
stop_condition=lambda y, t: y[0] < tphi - 2 * pi or
y[1] < 0)
return y[-1]
def find_gamma2(k0, k1est):
c.k[0] = k0
tphi = arccos(pd2.omega * pd2.omega / c.k[0])
z0 = integrate_for_gamma2(k1est)
k1start = k1est
k1end = k1est
if z0[0] > tphi - 2 * pi:
i = 1
while z0[0] > tphi - 2 * pi:
k1start = k1end - (powi ** i) * tol
z0 = integrate_for_gamma2(k1start)
i += 1
if i == imax:
print('Error: find_gamma2, imax too low!')
print('k0 = {0}, k1est = {1}'.format(k0, k1est))
return None
else:
i = 1
while z0[0] <= tphi - 2 * pi:
k1end = k1start + (powi ** i) * tol
z0 = integrate_for_gamma2(k1end)
i += 1
if i == imax:
print('Error: find_gamma2, imax too low!')
print('k0 = {0}, k1est = {1}'.format(k0, k1est))
return None
while k1end - k1start > tol:
z0 = integrate_for_gamma2(0.5 * (k1start + k1end))
if z0[0] > tphi - 2 * pi:
k1end = 0.5 * (k1start + k1end)
else:
k1start = 0.5 * (k1start + k1end)
return 0.5 * (k1start + k1end)
# Вычисление бифуркационных значений
k0 = linspace(pd2.omega * pd2.omega + 1e-4, 3, 100)
gamma1 = empty_like(k0)
gamma2 = empty_like(k0)
gamma1est = 0.256
gamma2est = 0.260
for i in range(k0.size):
gamma1est = find_gamma1(k0[i], gamma1est)
gamma1[i] = gamma1est
gamma2est = find_gamma2(k0[i], gamma2est)
gamma2[i] = gamma2est
# Пределы по осям
k0lim = (1., 3.)
k1lim = (-0.5, 1.5)
nk0 = 201
nk1 = 201
k0g = convolve(linspace(pd2.omega ** 2, k0lim[1], nk0), [0.5, 0.5], 'valid')
k1g = convolve(linspace(-pd2.lambd, k1lim[1], nk1), [0.5, 0.5], 'valid')
K0, K1 = meshgrid(k0g, k1g)
def find_maximal_re(k0_0, k1_0):
return amax(real(roots([1. - pd2.alpha, pd2.lambd + k1_0,
-pd2.omega * pd2.omega + k0_0])))
vfind_maximal_re = vectorize(find_maximal_re)
maxrez = vfind_maximal_re(K0, K1)
im = imshow(-maxrez, interpolation='bilinear', origin='lower', cmap=cm.gray_r,
extent=(pd2.omega ** 2, k0lim[1], -pd2.lambd, k1lim[1]))
im_bar = colorbar(im, orientation='vertical', shrink=0.8,
label=r'$-\max \,\mathrm{Re}\, p$');
plot([pd2.omega ** 2, pd2.omega ** 2, k0lim[1]],
[k1lim[1], -pd2.lambd, -pd2.lambd], 'pink')
plot([pd2.omega ** 2, k0lim[1]], [pd2.lambd, pd2.lambd])
plot(k0, gamma1, k0, gamma2)
axis([k0lim[0], k0lim[1], k1lim[0], k1lim[1]])
xlabel(r'$k_0$')
ylabel(r'$k_1$')
title(r'Бифуркационная карта');
savefig('files/bifurcation_map.png', dpi='figure', format='png',
transparent=False, frameon=False)
savefig('files/bifurcation_map.pdf', dpi='figure', format='pdf',
transparent=False, frameon=False)
